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ABSTRACT 

Recent numerical magnetohydrodynamic calculations by Braithwaite and collabora- 
tors support the 'fossil field' hypothesis regarding the origin of magnetic fields in 
compact stars and suggest that the resistive evolution of the fossil field can explain 
the reorganisation and decay of magnetar magnetic fields. Here, these findings are 
modelled analytically by allowing the stellar magnetic field to relax through a qua- 
sistatic sequence of nonaxisymmetric, force-free states, by analogy with spheromak 
relaxation experiments, starting from a random field. Under the hypothesis that the 
force- free modes approach energy equipartition in the absence of resistivity, the output 
of the numerical calculations is semiquantitatively recovered: the field settles down to 
a linked poloidal-toroidal configuration, which inflates and becomes more toroidal as 
time passes. A qualitatively similar (but not identical) end state is reached if the mag- 
netic field evolves by exchanging helicity between small and large scales according to 
an a-dynamo-like, mean-field mechanism, arising from the fluctuating electromotive 
force produced by the initial random field. The impossibility of matching a force-free 
internal field to a potential exterior field is discussed in the magnetar context. 
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1 INTRODUCTION 



Magnetars are neutron stars whose (un)pulsed X-ray luminosities ( < IQ-^^ W) exceed their s pin-down luminosit ies (IThompson fc Duncan 



19951 ). They fall into t wo classes: soft gamma repeaters (SGRs) (^Rothsc hild et al 



1994 



Hurley et al.l 2005) and anomalou s 
20031 : IWoods fc ThompsonI bood ). 



X-ray pulsars (AXPs) ([Thompson fc Duncan! Il996l : iMereghettil I200Q : .McLaughlin et al. 
From their observed spin periods and spin-down rates, one can infer that magnet ars possess very high surface magnetic 
fields B ^ 5 X 10^ 1 0, high er than most (although not all) rotation-powered pulsars ( Kouveliotou et al. 19981 : MelatosI 19991 : 
ji3l2 



Kaspi fc McLaughlinll2005l ). Duncan and Thompson (1992) proposed that the resistive decay of the magnetic field powers the 
continuous X-ray emission from AXPs and the intense outbursts (lasting a fraction of a second) fro m SGRs, e.g. Palmer et 
al. (2005). AXPs may be SGRs which are yet to burst ([Thompson fc Duncanll 19961 : iMereghettilboOOl ). 

It rema ins unclear how such strong stellar magnetic fields come to exist. Broadly speaking, there are two hypotheses for 
their origin ( Braithwaite fc Spruit|[2006l ) : (1) fossil field, where the progenitor star's magnetic field is amplified by compression 
during core collapse and is frozen into the highly conducting compact remnant; and (2) dynamo amplification, where the field is 
generated by a convective dynamo in the proto-neutron star. (In this paper, we discuss dynamos operating in the protostar stage 
rather than continuously regenerating the magnetic field of a neutron star, white dwarf, or Ap star during the star's lifetime; the 
latter idea seems to be ruled out by observations of Ap stars and white dwarfs, as we discuss below.) Observations to date do not 
favour conclusively one hypothesis over the other. Indeed, the hypotheses are not mutually exclusive; a dynamo operating in the 
protoneutron star/white dwarf stage may amplify the fossil field inherited from the main sequence. Fossil fields are supported by 



* E-mail: amastran@physics.unimelb.edu.au 
t E-mail: amelatos@physics.unimelb.edu.au 
^ Throughout this paper, we use SI units. 1 T 



= 10" G. 



2 A. Mastrano and A. Melatos 



the observed magnetic field distribution in Ap stars and white dwarfs ( Ferrario fc Wickramasinghel 200^ ). On the other hand, 
the magnetar birth rate implied by observations agrees with the operation of a large-scale dynamo in rapidly rotating, massive 
proto-neutron stars (with masses > 25M^ and periods ^ 1 ms) jKouveliotou et al. 1994 : Gaensler. Gotthelf. fc Vasisht 19991 : 
Heeer et alJbooi iGaensler et al.ll2005l : IPopov fc Prokhorovll2006l l. On the theoretical front, one of the major issues with the 



fossil field hypothesis is the (in)stability of the progenitor field. Much work has been done to show that purely poloidal and 
purely toroidal field configurations are unstable (Prendergast 1956; Wright 1973; Goossens & Tayler 1980; van Assche, Tayler, 
& Goossens 1982), but a linked toroidal-poloidal configuration may be stable (Tayler 1980). 

There are, in fact, two different fossil field hypotheses: "strong", where the magnetic field of the progenitor is ~ 1 T for a 
IMq star, following compression of a strongly magnetised (3 x 10~^ T) molecular cloud (Moss 2001, 2003), and "weak", where 
the field is generated after the progenitor forms by various dynamo processes [see e.g., Charbonneau fc MacGregor (1996) and 
Spruit (2002)]. In this paper, we do not distinguish between the two hypotheses, as our analysis is too elementary to treat 
the magnetization of the progenitor. Magnetic flux conservation implies that a magnetar progenitor needs to have magnetic 
fleld strength i? > 1 T (Woltjer 1964; Vink fc Kuiper 2006). The core and envelope of the progenitor are strongly coupled, 
angular momentum is rapidly transferred away by a stellar wind, and the core rotates rather slowly, for example P ~ 10 ms 
for 15 Mq stars and P ~ 4 ms for 35 Mq stars (Spruit 2002; Heger, Woosley, & Spruit 2005). This is slower than the P ~ 1-3 
ms needed for an efficient dynamo to operate (Duncan fc Thompson 1992; Thompson & Duncan 1995, 1996; Heger et al. 
2005; Vink 2008). Heger et al. (2005) calculated the efi'ects of magnetic braking on newly born neutron stars and concluded 
that P < 3 ms is achieved only in stars with M > 25Mq (Gaensler et al. 2005; Heger et al. 2005). The requirement of rapid 
rotation suggests a stringent absolute test of the dynamo hypothesis: are enough neutron stars born with such rapid rotations 
to explain the observed magnetar birth rate? Taking 25M0 to be the mass cutoff, Gaensler et al. (2005) estimated the birth 
rate of magnetars to be ^ 10% that of radio pulsars, consistent with current populations of AXPs and SGRs (Kouvelioutou 
et al. 1994; Gaensler et al. 1999). In contrast, the explosion energies of the supernova remnants associated with the AXPs IE 
1841-045 (Kes 73) and 1E2259+586 (CTB109) and the SGR 0526-66 (estimated from the shock speed, which is determined 
from X-ray spectra measured by XMM-Newton) suggest an initial spin period P > 5 ms, not enough for a dynamo (Vink & 
Kuiper 2006; Vink 2008). 

The lack of correlation between the observed field strength and spin periods of Ap stars and white dwarfs seems to argue 
against the hypothesis that a dynamo can operate in the core to regenerate the magnetic field continuously (Braithwaite 
& Spruit 2004; Petit et al. 2007). Observations of main sequence stars that are potential magnetar progenitors do not 
conclusively favour either fossil fields or protostar dynamos. From the similar fluxes of magnetic A stars (called Ap stars) and 
white dwarfs, Ferrario and Wickramasinghe (2006) inferred that O and B progenitors need to have dipolar fields of strength 
10"" T < B < 10"^ T to generate neutron stars with lO'' T < B < lO" T. Unfortunately, the magnetic fields of these O 
and B stars are hard to observe. This is because their absorption lines in the optical range are few and broad (perhaps due 
to rotation), making the Zeeman effect, by which the magnetic fields are measured, harder to detect (Donati et al. 2006a, 
b; Vink 2008). Two magnetic stars, HD 191612 (O class) and r Sco (B class), are observed to have fairly long periods, 41 
and 538 days respectively, which are inconsistent with the dynamo hypothesis (Donati et al. 2006a, b). Note, however, that 
their current rotation periods may not refiect directly on their periods during the protostar phase [many processes, such as 
accretion disk locking, can slow down the star after the stable field is in place (Braithwaite 2008)]. So, while the magnetic 
field today may not be generated continuously by a dynamo, it may arise from a dynamo during the protostar stage. It must 
also be kept in mind that our knowledge of magnetism in main sequence stars of masses > 8Mq is sketchy. 

Recent magnetohydro dynamic (MHD) simulations have shown that an unstable, rando m initial field can decay into a long- 
term stable configuration I Braithwaite fc Nordlundl 20061 : Braithwaite fc Spruitll2004l . l2006h . Starting from a random magnetic 
field, the simulations follow the field as it relaxes towards a lower energy state, with and without diffusion, transforming 
spontaneously into a stable configuration: a nearly axisymmetric torus with roughly equal poloidal and toroidal fluxes, which 
expands as a function of time. P oloidal field lines are wrapped around the torus, extending through the atmosphere and 
creating an 'approximate dipole' ([Braithwaite fc Nordlundl bood l . It must be noted that the results of the cited numerical 
simulations do not say anything directly about whether a neutron star's field comes from a dynamo process operating during 
the protoneutron star phase or simply conserved from the main sequence progenitor. A stable magnetic field configuration is 
necessary in both fossil field and dynamo hypotheses. 

Motivated by these numerical simulations, we ask here whether it is possible to model the magnetic evolution ap- 
proximately yet analytically as a slow, helicity-conserving relaxation through a quasistatic sequence of force-free states, by 
analogy with spherom ak experiments. A previous paper on relaxing force-free magnetic fields focused on axisymmetric modes 
i Broderick fc Naravani.2007 ). In this work, we draw particular attention to the importance of nonaxis ymmetric modes, without 



which it is impossible to reproduce the linked poloidal-toroidal field lines observed in the simulations IjBraithwaite fc Nordlund 



20061 : [Braithwaite fc Spruitl[2004l . 12006 [). We also allow our model neutron star to develop a small-scale turbulent velocity field 



as the magnetic fiel d relaxes, generalizing Br oderick and Narayan's (2007) work to accommodate (some of) the physics of a 
mean-field dynamo ( Blackman fc Fieldl 2004 ) . 

In Section 2, we delineate our aim and the four hypotheses regarding the nondiffusive evolution of the force-free field 
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that we investigate. In Section 3, we discuss the basic properties of force-free magnetic fields in detail. Then, in Section 4, 
we discuss the results of our study of the nondiffusive evolution of a composite force-free field through the four different 
hypotheses. In Section 5, we briefly discuss the diffusive evolution of a force-free field. Finally, in Section 6, we present our 
conclusions and discuss in greater detail the fundamental problem of matching a force-free magnetic field to a potential field. 



2 EVOLUTIONARY HYPOTHESES 

The aim of this paper is to construct a simple, approximate, analytic model for the evolution of magnetar fields (via non- 
diffusive and diffusive processes) , which reproduces coarsely the numerical results of Braithwaite and collaborators. In pursuing 
this objective, we hypothesize that the field evolves according to two principles: (1) it passes quasistatically through a series 
of force-free equilibria, each a superposition of axisymmetric and nonaxisymmetric modes; and (2) the relative weightings 
of the modes are controlled by the conservation (or slow evolution) of global quantities like magnetic helicity, energy, or 
hydrodynamic vorticity. Obviously, this approach sacrifices realism for simplicity; it is emphatically not a substitute for 
high-resolution, resistive MHD simulations. 

We choose force- free fields as our basis because of their astrophysical prevalence (Priest 19841 ) and analytic simpli 



ic- 



ity. This choice is justified by laboratory experiments with spheromaks, which are topologically congruent with a star, 
where it is observed that a spherical magnetic field relaxes into a force- free axisymrnetric state, known as the Taylor state 



( Rosenbluth fc Bussac 19791 : Goldenbaum et al. 1980l : Yamada et al. 1981 ; Janos et al.ll 19851 : lTavlor|[l986h . In a magnetar, the 



force-free approximation is defensible with respect to latitudinal and azimuthal force components, with Vp ^ /iQ"^(V x B) x B 
in these directions (where p is the pressure and B is the magnetic field). The approximation breaks down in the radial direction 
in a realistic star, due to the strong gravitational force. However, the main aim of this paper is to make contact with the 
MHD simulations of Braithwaite and collaborators. For this purpose, as noted below, the force-free approximation is adequate 
overall. 

We stress here that nothing physical requires the field to be force-free; indeed, with VP 2> J x B, one would naively 
expect B to be anything but force-free (at least with regard to the radial gradients). We want to make it clear that we do 
not propose force-free fields as a realistic model, but rather as a convenient mathematical artifice for studying the long-term 
evolution of magnetic helicity and field geometries in a neutron star. That is, we take the simplifying assumption that the 
field is force-free, to see whether the behaviour of such a field mimics that of a more general field investigated by Braithwaite 
et al. Additionally: 

(i) Force-free fields are an analytically convenient orthonormal basis in which any interior field can be expanded. It will 
be worth testing empirically (from the simulation output) what the coefficients in a force-free expansion are, in view of the 
qualitative visual agreement with the simulation results reported in Sections 4 and 5. 

(ii) In many other high-/3 (|Vp| ^ |J x B|) plasmas, e.g., in spheromaks (Jarboe et al. 1980; Turner et al. 1983; Janos et 
al. 1985; Barnes et al. 1986; Knox et al. 1986), the field does relax to a force-free (minimum energy) state (Taylor 1986), even 
though it does not have to and one naively would not expect it to. It does this by tending "thermodynamically" to a state 
with Vp ~ and J x B ~ separately (because a magnetic field does no work in the first law of thermodynamics). The 
analogue here is a state with — Vp + pV"!? ~ and J x B ~ 0. 

The MHD simulations cited above start from a random field. The random field is not force free, but this is immaterial: the 
role played by the initial state, in our model, is to define the initial values of conserved (or nearly conserved) global quantities 
like magnetic helicity. We focus on the three lowest poloidal orders of the magnetic field, including nonaxisymmetric modes. 
Our initial state is a sum of these modes, equally weighted to emulate a random field (i.e., spatially uncorrelated) where all 
modes (up to some wavenumber cutoff) are equally probable. 

Next, we entertain four possibilities for how the field evolves without diffusion, as in the first stage of the numerical 
simulations of Braithwaite and Nordlund (2006). One possible end state of nondiffusive evolution is energy equipartition, 
where each force-free mode contributes equal energy to the total. Sec. 4.2 explains the equipartition assumption in more 
detail. We conjecture that a random field may evolve into an equipartitioned state either by keeping the total energy and 
helicity constant (state A) or by keeping the total helicity constant while allowing the energy to vary (state B). A second 
possible end state of nondiffusive evolution is dictated by the dynamo theory of Blackman and Field (2004), in which helicity 
is exchanged between short- and long- wavelength modes through the agency of a turbulent electromotive force (a effect). 
Again, the end state is reached either by keeping the total energy and helicity constant (state C) or by keeping the total 
helicity constant while allowing the energy to vary (state D). The discretisation scheme of the simulation gives rise to a small 
"effective diffusivity", which leads to some dissipati on and energy loss (by a factor of ~ 10 - 100) during the initial relaxation, 
before a stable magnetic field configuration forms ( Braithwaitg 12008 ) . In our analysis, however, we prefer to keep an open 



mind and allow cases without (A and C) and with (B and D) energy loss. 

Lastly, we briefiy discuss how the relaxed, force-free field diffuses over long times when the electrical resistivity is not 
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zero. In the numerical simulatio ns, the magnetic field is observed to diffuse out of the star in a roughly dipolar configuration 
( Braithwaite fc Nordlundlbood ). We show that a field which is a linear combination of several force-free modes diffuses while 
keeping the force- free parameter a (defined in Section 3) and relative weighting of the modes constant, thus retaining its 
pre-diffusion shape. We also discuss the difficulty of matching a force-free field with spatially uniform a to a potential field 
outside the star. 



3 NONAXISYMMETRIC FORCE-FREE FIELDS 



Force-free magnetic fields obey the condition V x B = oB, where a(x) (units: inverse length) is a function of position in 
general. The conduction current density, J — fJ-Q^V x B, is therefore parallel to B and the Lorentz force, proportional to 
J X B, vanishes. Hence, force-free fields offer an attractive, static, long- term magnetic configuration. In generah the solution 
to the force-free condition that also fulfills V ■ B = can be written as (jChandrasekhar fc Kendallll 19571 : |Priestlll984r ) 



B = aV X (Tpn) + V X [V X (Tpn)] (1) 
where n is some vector and t/j is a scalar solution of the Helmholtz equation 

(V^ + a^)i^ = 0. (2) 

choosing n = and restricting attention (for simplicity) to uniform a, we can 



In spherical polar coor dinates (r, 9 , 
write the magnetic field as l|Priestlll984l ) 



Br = y^J2 Cnmn{n+l)r ^/2j„+i/2(ar)Pr(cos6i)exp(im(; 



n— 1 m— 



(3) 



Bit, 



-3/2 r , ,,dP^(cose) 



ji — 1 m— 



-1/2 



Jn+i/2{ar)Pn{cose) exp (im4>)/ f 



mtmmr 



-3/2 T c ^ 



-1/2 



Jn-i/2{<i''')]P^ {cos 9) exp (jm(/>)/ sin ( 



Cn 



-1/2 



+i/2{ar\ 



rfPr(cos6>) 
dB 



exp {imij. 



(4) 



(5) 



where Cnm are complex constants, J„±i/2 denotes Bessel functions of order n ± 1/2, and denotes associated Legendre 
functions of the first kind. Physical field components are the real parts of (3)-(5). We plot some of these modes in Fig. 1. The 
leftmost plot in the top row shows the axisymmetric Taylor state, the end state achieved by relaxing spheromak magnetic 
fields (Taylor 1986). Note that a nonaxisymmetric vector field may have chaotic (albeit closed) field lines (e.g. second row, 
centre plot, Fig. 1) e ven though the vector field is smooth, because the field line equations form a nonlinear dynamical system 
of dimension three ( Botha fc EvangelidisI 2004 ). The field lines are closed (because V ■ B = 0) but chaotic (two field lines 
that are arbitrarily close at some point will diverge macroscopically along their length). By contrast, the field lines of an 
axisymmetric vector fie ld must form simple loops, by the Poincare-Bendixson theorem for two-dimensional dynamical systems 
( Jordan fc Smithll 19991 ). 

A nonuniform a gives an extra term in the diffusion equation, since j^V^B — —a^rjQ — i]{Va) x B. Under the action of 
the new term, the field changes shape as it diffuses. In addition, a nonuniform a facilitates global current closure, since the 
currents carried by two flux tubes of equal area (proportional to the local value of a) are not necessarily equal if a varies 
from place to place. This allows return currents to be confined to thin sheets that leave most of the field undisturbed, helping 
greatly with the matching problem discussed in Sec. 6. 



3.1 Helicity 

The magnetic helicity within a volume V is defined as 



H : 



dVA ■ B, 



(6) 



where A is the vector potential. It has been argued th e oretically (Tavio3 19861) . and observ ed in spheromak experiments 



i Jarboe et al]|l980l : [Turner et aLlll983l : iJanos et al.lll985l: iBarnes et al.l 



1986 



Knox et al 



19861 ). 



a force- free field relaxes. More generally, if a nonzero electric field E 



-V0 



that H remains constant as 
dA/dt is included (where is the scalar 



(1,0) 




Figure 1. Magnetic field lines for the force-free modes l^n^3, O^m^n, with the magnetic foot points (the starting points for the 
field line integration) initiahsed at ars = 4.49, arccos(0) = (—8/9, —7/9, 7/9, 8/9), (j> = (0, 7r/8, it/A, 157r/8). Each panel is labelled 
by (n,m). The north-south axis lies vertical in the plane of the page. 



potential), Faraday's law and the definitions of H and A imply dH/dt = d(A ■ B)/dt = 
one can then integrate over a volume V to obtain (jBrandenburg fc Subramanianll2005l ) 



dH 



= -2 



dVE ■ B 4- / dS ■ (A X E + 0B). 



av 



-2E ■ B + V ■ (0B + A X E), which 



(7) 



The helicity is gauge-independent if B is strictly tangential on the bounding surface dV , which is true for spheromaks but not 
for a neutron star, unless we extend dV to infinity. This makes sense only if B falls off fast enough, such that H is finite. To 
achieve this, the force-free field [which scales as |B| oc for large r, according to (3)-(5)] must be matched to a potential 
field at r = (where is the surface of the star), something that is difficult to do both in practice and in principle (Sec. 
6 further explains why) ( Wheatlandl 2006[ ). Consequently, in this paper, we assume that the field is force-free everywhere 
while integrating over the volume r ^ in (6) to keep H finite. The results in Sections 4-5 must be interpreted with this 
approximation in mind. 

The reader should note that, for simplicity's sake, we implicitly assume that the field at r > is also described by Eqs. 
(3)-(5), not by a potential field. This means that, at infinity, the magnetic field does not approach a dipole field (which is 
physically unrealistic and contrary to pulsar observations). Ideally, we wish to use an external field which is more realistic, 
such as a dipole field or a dipole-dominated multipole field, but this raises the issue of boundary matching, which is nontrivial 
and is discussed further in Sec. 6. 



3.2 Energy 

In a force-free field with uniform a, we have H = Jy dVA ■ B — (l/a) J-^ dVB^ — 2iioW/a, where W is the total magnetic 
energy and is the magnetic permeability of free space. This is true both for each mode (n, m) in (3)-(5) and for a linear 
superposition of modes; it can be easily verified by direct integration that the cross-terms Jy dVJi(^i-^^rni) '^(i2,m2) with h ^ h 
and/or mi 7^ m2 all vanish because of the linear independence of the modes. For a particular value of a, higher mode numbers 
have higher energies. Fig. 2 illustrates this trend for modes with 1 ^ n ^ 3 by graphing W against a for fixed r^. The modes 
(n, m) — (1,0) and (1, 1) have equal energies, despite having different structures, because the energies in their radial and 
tangential components are separately equal. 




Figure 2. Magnetic energy Wnm for the force-free mode (n,m) in (3)-(5) (in units of C'^^/^q as a function of ars, where a is the 
force-free parameter and Ts is the radius of the star. SoUd curves: n = 1; dashed curves: n = 2; dashed-dotted curves: n = 3. For each 
value of n, the lowest curve corresponds to the axisymmetric (m = 0) and the highest curve corresponds to the m = n mode. Note that 
Wio = VKii and Wao = W31. 



Taking the dot product of Faraday's law with B/(2/io), then using E = J/g — v x B to ehminate E, and fina lly integrating 
over a volume V , one obtains the following equation for energy evolution ( Brandenburg &: Subramanian 20051 ) . 



dW 
dt 



[ dV{J X B) • V - i / dVJ^ - — [ 
Jv Jv MO Jdv 



rfS ■ E X B, 



(8) 



where v is the velocity of the charged fluid threaded by B, cr is the conductivity, and J is the current density. For force- 
free flelds, one can write fioJ = V x B = aB. Thus, the flrst term on the right-hand side of Eq. (8) vanishes, giving 
ExB = (J/cr — vxB)xB = — (vxB)xB and hence 



dW 

dt 



2 1 
B^dV+ — 



MO Jv 



dS ■ \(v x B) x Bl 



(9) 



4 PRE-DIFFUSION EVOLUTION 

The numerical experiments of Braithwaite and collaborators ( Braithwaite fc Nordlundl2006l : Braithwaite fc Spruitll2004l . l2006l ) 
proceeded in two stages. In the first stage, a random magnetic field threading a stationary ball of plasma is loaded into a 
three-dimensional MHD solver with both numerical and physical resistivities. We study the ideal-MHD evolution in this 
section. In the second stage, the collisional resistivity is switched on. We study the resistive evolution in Section 5. Note 
that we do not mea n to imply here that the only diffusion process is Spitzer resistivity [another possibility is, for example, 
ambipolar diffusion (|Goldreich fc ReiseneggeJl 19921 )]. Our results apply regardless of the actual diffusion mechanism, provided 
it is isotropic (i.e., cr is a scalar). 



4.1 Initial conditions: random field 

Braithwaite and Nordlund (2006) initialised their simulation with a random field, whose maximum wavenumber in a Fourier 
decomposition is set by the grid scale (~ lO^iir^)- We simulate the initial state as a linear combination of the modes 1 ^ n ^ 3, 
^ m ^ n, in which the modes are equally weighted, i.e. all constants Cnm are set to be equal. We draw the field lines of 
the initial state in Fig. 4a. Note that the maximum wavenumber is fcmax ~ 4.27r (for n = m — 3). One can readily increase 
fcmax by increasing n to match Braithwaite and Nordlund's (2006) experiment, but there is no need; the conserved global 
quantities are approximated to ~ 10 % with n ^ 3. 

Note that the total energy Wt = X]n=i X^m^o '^"m for our initial state is minimised at a = 7.0. 
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Figure 3. Hclicity [solid curve, in units of (Cio)^/a] and its derivative witli respect to a [daslied curve, in units of (Cio)^/a^] versus 
ars, for the composite equipartitioned state witli n = 1, 2, 3 and ^ m ^ n. 



4.2 Relaxation into equipartition 

We now test the hypothesis that the force-free magnetic field evolves to energy equipartition. That is, we suppose that the 
relative weightings (C„m) of the modes adjust until each mode contributes equally to the total energy, Wt = A''(A''-|-l)Wio (for 
1 ^ n ^ A?^, ^ m ^ n). In ideal MHD, this cannot happen resistively. We postulate instead that it happens in the simulations 
because the random initial perturbations 5B (and concomitant fluid velocity perturbations 5v) consp ire to produce a nonzero 
mean electromotive field (Jv x (5B) 7^ as they "unwind", just like the a effect in mean- field MHD (jBiskampllloii ). 



A force-free magneti c field relaxi ng via the a effect seeks out a state which minimizes its helicity at a nonzero value, as 
observed in spheromaks ( Taylor 19861 ). We find dH/da = at ar^ = 4.49 (where is the radius of the model neutron star), 



the same conclusion reached by Taylor (1986) (see Fig. 3). This result applies to any combination of modes in equipartition, 
because Wt is proportional to Wio, and the energy of the (1,0) mode is minimized at avs = 4.49. Hence, Taylor's minimum 
eigenvalue for axisymmetric spheromak eigenfunctions applies just as well to nonaxisymmetric eigenfunctions if the modes 
tend towards equipartition. We plot the helicity and its derivative for the composite equipartitioned state as functions of a in 
Fig. 3; the local minimum at a = 4.49r7^ is apparent. The helicity is calculated over the stellar volume, not out to infinity. 
This choice is bound up in the physics of the boundary conditions, which is discussed further at the end of this section and 
in Section 6. 

The total helicity Ht is related to the total energy by Ht = 2/ioWT/ffl. Therefore, if the field is to change from a random 
configuration (where all Cnm are equal) to an equipartitioned one (where all modes contribute equal energy) while keeping 
both helicity and energy constant, a must remain constant as well. We designate this final state as State A. If, on the other 
hand, the field is allowed to convert its energy into, say, kinetic energy (e.g., random fiuid motions in the MHD simulation), 
a can change to minimize Wt independently. We designate the corresponding final state as State B. We present the values of 
the mode coefficients C„m. achieved through these two scenarios in Table 1. Note that, for all four states, all constants Cnm 
increase from their initial values except for the (3,2) and (3,3) modes. Although Cnm takes higher values in State B than in 
State A, State B actually has a lower magnetic energy, due to the decrease in a. 



We plot magnetic field lines of the initial state. State A, and State B in F igs. 4a, 4b, and 4c respectively . The two final 
states are similar in shape, with State B having a larger radius tm, defined as ( Braithwaite fc NordlundluOOfil ) 

'^M = Vttt^ • (10) 



As the star does not have an actual surface, we use tm as a proxy for the characteristic spatial extent of the field. Specifically, 
we find rivi/j's = 0.71, 0.85, and 1.0 for the initial state. State A, and State B respectively. Fig. 4 is to be compared with the 
bottom left panel of Fig. 6 in Braithwaite and Nordlund (2006), which depicts the state achieved by the ideal-MHD simulation 
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Table 1. Parameters and constants of the initial and final states of the force- free magnetic field, under four evolutionary scenarios 
without diffusion. The force-free parameter a is measured in units of rj^. The constants Cnm are measured in units of Cio at the 
start of the pre-difi'usion phase, denoted by C. 



Scenario 


ClTs 


Cio/C 


Cii/C 


C20/C 


C21/C 


C22/C 


C30/C 


C31/C 


C32/C 


C33/C 


Initial state 


7.0 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


State A: equipartition. 


7.0 


11 


11 


8.1 


4.7 


2.3 


7.0 


7.0 


0.90 


0.37 


//t, Wt, and a constant 






















State B: equipartition, 


4.5 


11 


11 


8.4 


4.8 


2.4 


8.8 


8.8 


1.1 


0.47 


constant 






















State C; mean-field MHD, 


7.0 


20 


14 


8.0 


4.1 


1.6 


5.6 


5.3 


0.60 


0.21 


Ht, Wt, and a constant 






















State D: mean-field MHD 


4.5 


21 


15 


8.1 


4.2 


1.7 


5.7 


5.4 


0.61 


0.21 



i?x constant 



(a) 




Figure 4. Magnetic field configurations produced by the evolutionary hypotheses in Section 4.2. (a) Initial random field, with Cnm = 1 
for all 1 ^ n ^ 3, ^ m ^ n, and ar^ = 7.0. (b) Final force-free field in State A, with Cio = 11, Cn = 11, C20 = 8.1, C21 = 4.7, 
C22 = 2.3, C30 = 7.0, C31 = 7.0, C32 = 0.90, C33 = 0.37, and avs = 7.0. (c) Final force-free field in State B, with Cio = H, Cii = H, 
C20 = 8.4, C21 = 4.8, C22 = 2.4, C30 = 8.8, C31 = 8.8, C32 = 1.1, C33 = 0.47, and ar^ = 4.5. The shaded sphere is at a radius r^. 



(a) 
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(b) 





Figure 5. Magnetic field configurations produced by the evolutionary hypotheses in Section 4.3. (a) Final force-free field in State C, 
with Cio = 20, Cii = 14, C20 = 8.0, C21 = 4.1, C22 = 1-6, C30 = 5.6, C31 = 5.3, C32 = 0.60, C33 = 0.21, and ars = 7. (b) Final 
force-free field in State D, with Cio = 21, Cn = 15, C20 = 8.1, C21 = 4.2, C22 = 1-7, C30 = 5.7, C31 = 5.4, C32 = 0.61, C33 = 0.21, and 
avs = 4.5. The shaded sphere is at a radius r^. 



after 9 Alfven crossing times of nondiffusive evolution. Qualitatively, the results are similar. Note how field lines form toroidal 
loops around the star in both the numerical and analytical calculations. The toroidal loops are wound closer to the star in 
Fig. 4c than in Fig. 4b, indicating that State B resembles Braithwaite and Nordlund (2006) more closely than State A. Also 
note that the "surface" of the star is located at r = (where we also initialise the field lines), unlike in Fig. 6 of Braithwaite 
and Nordlund (2006), where the sphere is at r = O.Sra. Our model star does not have a hard boundary, and we do not match 

the "interior" force- free field to an "exterior" potential field at r — Vs- 

Imposing a boundary condition at r — rs, such as a conducting wall (more a ppropriate for laboratory plasmas') JSalingarosI 

199C ) , balance between external plasma pressure and interna l magnetic pressure IZhugzhdallQoil : Vandas. Romashets. fc Wataril 
2005 ). and force balance with an external magnetic field ( Low fc Bergen l2003l ) . affects the relaxation path and final state 



I Tavlon 1986: Botha fc Evangeli dis 2004). The problem of matching a force-free field to an external field is nontrivial 
Broderick fc NaravanI 2007 ). requiring a to vary (dis)continuously along the interface (plus the existence of non-force- free 
regions in the star). This matter is discussed in more detail in Section 6. Another potential issue is that, while the conductivity 
certainly damps the higher modes more quickly in a realistic field, as in the simulations, this is not true for the special case 
of force-free fields, where V^B = — a^B for all (/,m). This is another argument against the realism of force-free fields. 



4.3 Blackman and Field (2004) mean-field dynamo 

In the previous section, we do not advance a theory for how the parameter a and the constants Cnm change with time. We 
simply postulate an end state which fulfills energy equipartition and which is linked to the initial state by globally conserved 
quantities. In this section, we model the evolution more closely by applying a model developed by Blackman and Field 
(2004). This model is intended by the authors to be applicable to both astrophysical phenomena (such as stellar coronae) and 
laboratory plasma experiments. Blackman & Field (2004) explicitly sidestepped discussion of physical boundary conditions 
(particularly current closure) by utilizing periodic boundary conditions instead. 

The model, rooted in the principles of mean-field MHD, gives a recipe for the coupled evolution of the spatially-averaged 
helicity on two length scales (one large, one small), as well as the spatially-averaged fluid vorticity and electromotive field. 
Upon splitting the magnetic fleld, vector potential, and current density into (slowly varying) large-scale and (turbulent) 
small-scale components, the volume- averaged induction equation and Faraday's law imply 

dtHi^2{e-B)-2ri{B-J) (11) 
and 



9tJ/2 = 2{e-B)-2j7(b-j>. (12) 
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Table 2. Poloidal flux, toroidal flux, and ratio of poloidal flux to toroidal flux, for the random initial state and States A-D. The 
fluxes are expressed in units of CrJ^, where C is the initial value of the coefficient Cio- 



State 


Poloidal flux 


Toroidal flux 


Poloidal flux/toroidal flux 


Random state 


1.3 


2.6 


0.50 


State A 


14 


20 


0.70 


State B 


34 


31 


1.1 


State C 


27 


22 


1.2 


State D 


22 


47 


0.47 



In (12)-(13), rj is the magnetic diffusivity. Hi and H2 are the large- and small-scale magnetic helicities respectively, B and b 
are the large- and small-scale magnetic fields, J and j are the large- and small-scale current densities, and e — —E + rjCB-J) is 
the turbulent electromotive field. By utilising Ampere's Law, Faraday's Law, and the incompressible Navier-Stokes equation, 
Blackman and Field (2004) derived the following closed set of evolution equations for the spatially averaged helicities {Hi 
and H2), vorticity Hv = (v • V x v), electromotive field Q, and energy e: 

^ = ~2Qh\'\ki/k2f''' - 2hi{ki/k2f/RM, (13) 

dh2 
dr 



= 2Qh\'^{ki/k2f'^ - 2h2/RM, (14) 



^ = _2(l - ki/k2)(ki/k2)^^^Qhy^ - 2K/Rv, (15) 
^ = ~{ki/k2f'^h\^\l/2,){h2 - K) + {ki/k2f'^h\'\/i - e^'^fQ, (16) 
1^ = -2(1 - ki/k2){ki/k2f'^Qh\''' ~ 2e/Rv. (17) 



In (14)-(18), we define k\ and as the characteristic large- and small-scale wavenumbers, Rv = H2{0y'^'^ /vk\^'^ (where 
u is viscosity), Rm ~ H2{QY^^ /riklj"^ , = Hv /k2H2{0), hi and /12 are the large- and small-scale helicities normalised to 
1/2(0), Q — ~e\\/k2H2{0) is the component of the turbulent electromotive field parallel to B, e = (v^)/fe2//2(0), where v is 
the turbulent fiuid velocit y, t = tk^^^ H2(0Y^^ , an d / ~ 1 is a dimensionless constant which parametrizes the (unspecified) 



microphysical dissipation (|Blackman fc Fieldll2004l ). The notation '(0)' denotes an initial value 



Let us now look at what happens in the steady state. From (14) and (15), the steady-state helicities satisfy 

h2 = {ki/k2fhi. (18) 

For our force-free modes (3)-(5), we can write kmn oc (n^ + m^y^^. Suppose that the two-mode result (19) can be generalized 
to apply to all the modes in a linear combination. We designate the final relaxed state with C„m related through (19) as 
State C if both helicity and magnetic energy are kept constant throughout relaxation, and as State D if only the helicity is 
kept constant. Magnetic field lines for States C and D are shown in Figs. 5a and 5b respectively. These states do not obey 
the equipartition relation. Note that States C and D are very similar (see also Table 1), while A and B are less so, indicating 
a weaker dependence of Ht on a. The values of ryi/ra are 0.77 and 0.78 for States C and D respectively. As evident from 
Table 1, all Cnm in States C and D are smaller than in States A and B, except for Cim. In fact, the axisymmetric (1,0) 
mode dominates in States C and D, being approximately 20 times greater than in the initial state, giving these states a more 
toroidal shape (see Figs. 5a and 5b). 

To quantify further the 'shape' of the field, we compute the poloidal flux $e and toroidal flux <I>0 from 

$8=/ d(j> dr rBeir,n/2,<j,) (19) 
Jo Jo 

<i>^ = - de dr rB^{r,e,0) + - dO dr rB^{r,e,-R). (20) 

^ Jo Jo ^ Jo Jo 

The results are listed in Table 2. In all four final states, <^e and $0 are comparable. This is a natural property of the linked 
toroidal-poloidal field structure necessary for long-term MHD stability. Again, we see that nonaxisymmetric force-free modes 
play a key role, this time in ensuring long-term stability. 
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Figure 6. Dimensionless helicities of the (1,1) mode {hi, thick soUd curve), (3,3) mode (^2, thin solid curve), and vorticity {h^, dashed 
curve) as functions of r. h\ and h2 are initially equal (the modes are in equipartition) and hy is initially zero. The dimensionless physical 
parameters /, -Rjv/, and Ry are taken to be 1, 10, and 10 respectively. 



5 DIFFUSIVE EVOLUTION 

In the second stage of the numerical experiments described by Braithwaite and Nordlund (2006), the resistivity is switched 
on and the magnetic field obtained at the end of the (ideal-MHD) first stage is allowed to relax by Ohmic diffusion 



( Braithwaite fc Spruit 20061 ) . What are the implications for our analytic model? 



If the field remains in a linear superposition of quasistatic force- free modes, the linear independence of the modes ensures 
that they all decay at the same rate, proportional to the electrical conductivity a: 

dCnm/dt = {a^ / ^oa)Cnm- (21) 

If the force- free parameter a also remains constant throughout the diffusion process, it is evident that the structure of field lines, 
e.g. the quantity tm and poloidal-to-toroidal flux ratio, also remain unchanged. This is contrary to the results of Braithwaite 



and Nordlund, who found that the field becomes more poloidal and dipole-like as it diffuses outwards IjBraithwaite fc Nordlund 



200d ). By implication, the field does not remain force free (with constant a) in the resistive stage of the simulations. We note 



that an extra term is added to the magnetic field diffusion equation if a is nonuniform, since the term rfsJ^Ti now becomes 

— a^r;B — r?(Va) x B, which affects the shape of the field as it diffuses, depending on Vo; 

However, if we allow the energy to vary freely and even allow some loss of helicity ( Braithwaite fc Spruit 20061 ) ■ some of 



the magnetic energy at the start of the resistive stage is converted into fluid kinetic energy e, inducing a turbulent dynamo 
which regenerates the field (or, at least, exchanges helicity between small and large scales). Following Blackman and Field 
(2004), we consider two force-free modes with sufficiently different characteristic wavenumber k [modes (1,1) and (3,1), for 
example] and let them evolve from equipartition as per (14)-(18). The results are plotted in Fig. 6 and reveal some general 
trends. First of all, we find that the turbulent vorticity hv, which starts at zero, increases briefly then dies away, controlled 
by Rv (hv falls away slower for larger Ry, i.e. smaller viscosity). Secondly, the long- wavelength mode increases at first but 
then decreases while the short- wavelength mode decreases immediately (and faster). The rates at which the modes decay are 
controlled by the quantity Rm\ the larger the magnetic diffusivity, the faster the decay. Thus, in the long term, the lower 
n and m modes dominate, until all the higher order and nonaxisymmetric modes vanish; a force-free field which is a linear 
combination of several modes can decay into an axisymmetric Taylor state, given enough time and a mediating turbulent 
vorticity. We stress again, however, that this dynamo picture is one of many possibilities, and it does not take into account 
superfiuid and stratification effects in a neutron star. 



6 DISCUSSION 

We show in this paper that an initial magnetic field which is a linear combination of equally weighted, force-free modes can 
relax without diffusion, preserving magnetic helicity, to an equipartitioned state which possesses a linked poloidal-toroidal 
geometry (with toroidal flux dominating), just like the end result of the first (nondiffusive) stage of the numerical simulations by 
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Braithwaite and collaborato rs (Braithwaite fc Nordlund 20061 : Braithwaite fc Spruit 20061 ) ■ Alternatively, a turbulent, mean- 



field MHD dynamo model (jBlackman fc Fieldll20oi ) can be used to prescribe the evolution of the force-free field, allowing 



helicity to be exchanged between small and large scales through the agency of a turbulent electromotive field (and growth of 
small-scale hydrodynamic vorticity). Applying this particular model to our initial state, we again obtain a final configuration 
which qualitatively resembles the numerical results of Braithwaite and Nordlund (2006), but with smaller root-mean-square 
size rm. The inclusion of nonaxisymmetric force-free modes is important to obtain good agreement. 

A force-free field decays uniformly when a uniform resistivity is switched on, keeping the relative weighting of its modes 
(and hence its three-dimensional structure) constant. This beh aviour is not observed in the numerical simulations, wher e 
instead the field diffuses outwards over time into a dipolar shape ( Braithwaite fc Nordlundl 20061 : Braithwaite fc Spruit 20061 ) . 
If, however, a turbulent electromotive force is allowed to operate, we show that the higher-order modes decay faster, and the 
magnetic field tends overall towards the Taylor axisymmetric state (which is dominated by toroidal fiux at avs = 4.49). 

We conclude by speaking briefly about the vital issue of boundary conditions at the stellar surface. In this paper, we 
assume (for the sake of analytic simplicity) that the magnetic field is force-free everywhere, because our main goal is to 
elucidate the central role played by the nonaxisymmetric modes. In order to preserve finite helicity and energy, we are forced 
to integrate (artificially) A • B and /2^q over just the volume of the star, rather than over all space. Hence the helicity 
and energy we compute are semiquantitative approximations. It is natural to try and correct this shortcoming by matching 
a force-free field in the stella r interior to a current-free (potential) field in an insulating atmosphere above the stellar surface 
( Broderick fc NaravanlboO?! ). The matching conditions are popularly taken as: (i) normal component of B is continuous at the 
surface (which uniquely determines the external field) and (ii) the jump in the tangential component of B is proportional to 
the surface (skin) current (and indeed uniquely determines the latter). However, this prescription hides a subtle flaw: current 
continuity (V • J = 0) is violated at the surface, because it is impossible, in general, for the skin current leaving (entering) 
a surface patch sideways to balance the volume current fj,Q^aBr entering (leaving) the patch from below. The root of this 
apparent paradox concerning current closure can be traced back to the following classic theorem of magnetostatics: a magnetic 
field which i s force free everywhere within a finite volume cannot match to an external field which is potential everywhere 
I Priest 1984| 'FI. [P hysically, solutions with surface currents are not appealing, because conductivities drop off towards the 
surface of the star ( Braithwaite 20081 ) . In this way, Broderick and Narayan's (2007) treatment is problematic; surface currents 
are an artifact arising out of the matching problem and are accompanied by an artificial surface J x B force.] . One way to resolve 
the paradox is to allow currents to flow in the atmosphere above the surface, e.g. through discrete structures like magnetic 
arcades (cf. the Sun). Under these circumstances, the force-free parameter a(x) must vary with position, and the matching 
conditions at any interface r educe to continui ty of (i) the normal component of B(x), and (ii) a(x) (o ver one polarity of B), 
which are solved iteratively ( Wheatland 20061 ) . Loss of helicity through parts of the surface can occur ( Braithwaite fc Spruit 
20061 ) ■ although care must be taken to use appropriate boundary conditions on the slow relaxation time-scale. Another way 



to resolve the paradox is to create a non-force-free (i.e., stressed) region within the star, through which the currents can close. 
Imposing boundary conditions at the stellar surface which satisfy V ■ J = in this scenario is a nontrivial task. 

Both resolutions of the matching paradox outlined above may carry important implications for magnetar astrophysics. 
The first suggests the likelihood of magnetic activity (e.g. reconnection between magnetic arcades) in the magnetosphere, 
accompanied perhaps by radio and X-ray bursts. The second suggests that magnetic stresses may cause crust cracking in 
non-force-free regions. Note that neutron star magnetosphere models have popularly included surface currents in the past 
[e.g., Goldreich and Julian (1969), Ruderman (1975), Mestel et al. (1985), and Arons (1993)]. This paper exposes some of 
the limitations of analytic methods. A better way to study such possibilities self- consistently is to extend the simulations in 
Braithwaite and Nordlund (2006) and elsewhere to treat in more detail the mechanical behaviour of the stellar interior. 
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The theorem follows from the identity 

f dVx-(J X B) = — f dV\B\^ + — f dS- [2(x ■ B)B - IBI^x], (22) 
J J 2/10 J 

if we consider a spherical surface of radius R that fully encloses the force-free region and let /? — > oo. 
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